Baby Mental Life: Study 3

Preregistered analyses

# load required libraries
library(tidyverse)

── Attaching packages ────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.1     ✔ purrr   0.2.5
✔ tibble  2.1.1     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0
package ‘ggplot2’ was built under R version 3.5.2package ‘tibble’ was built under R version 3.5.2── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(langcog) # source: https://github.com/langcog/langcog-package

Attaching package: ‘langcog’

The following object is masked from ‘package:base’:

    scale

library(psych)

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

library(lme4)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following object is masked from ‘package:tidyr’:

    expand

library(kableExtra)
# set theme for ggplots
theme_set(theme_bw())
chosen_rot <- "oblimin"

# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/data_prep.R")

NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"
Joining, by = "question"
Joining, by = c("ResponseId", "attn_free_coded")
Column `ResponseId` joining character vector and factor, coercing into character vectorColumn `attn_free_coded` joining character vector and factor, coercing into character vectorJoining, by = "ResponseId"
Joining, by = "ResponseId"
NAs introduced by coercionJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vectorJoining, by = "capacity"
Column `capacity` joining character vector and factor, coercing into character vector

“Baby Mental Life: Study 3” was conducted on MTurk on 2019-04-24.

Our planned sample was 300 participants; based on Study 2, we initially recruited 352 participants. After filtering out participants who failed at least one of our attention checks, we ended up retaining 279 participants (retention rate: 79.3%), but at this point we have not supplemented with additional recruiting. At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.

In the end, we ended up with a sample of 279 participants who passed our attention checks, 249 of whom came from unique GPS coordinates (89.2%).

For this first pass, these data include participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.

Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 8 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable). In contrast to Study 2, participants completed all 13 ratings of a particular mental capacity on a single “trial” (rather than completing all 8 ratings of a particular target age on a single trial). In addition, participants answered questions about the “developmental factors” that might contribute to development in this domain.

For more details about the study, see our preregistration here.

Study 1 EFA

# load in EFA results from study 1
efa_S1 <- readRDS("../study 1/s1_efa.rds")

heatmap_fun(efa_S1) + 
  labs(title = paste0("STUDY 1 Parallel Analysis (rotation: ", chosen_rot, ")"),
       subtitle = "'% var.' indicates the amount of shared variance explained (total = 100%)")

Joining, by = "capacity"
Joining, by = "factor"

Plots: Capacity ratings as a function of target age and domain

ggplot(d_cap_rating %>%
         arrange(ResponseId, domain, capacity, target_ord) %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_),
           capacity = gsub("_", " ", capacity)),
       aes(x = target_num, y = response, group = ResponseId,
           color = domain)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  geom_path(alpha = 0.1) +
  geom_smooth(aes(group = capacity), 
              method = "loess", 
              color = "black") +
  scale_x_continuous(breaks = as.numeric(levels(factor(d_cap_rating$target_num))),
                     labels = levels(d_cap_rating$target_ord)) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Perceptions of development, by domain",
       subtitle = "Black line = loess smoothing",
       x = "Target age (numeric)",
       y = "Response (0 = not at all capable, 100 = completely capable)")

ggplot(d_cap_rating %>%
         arrange(ResponseId, domain, capacity, target_ord) %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_),
           capacity = gsub("_", " ", capacity)),
       aes(x = sqrt(target_num), y = response, group = ResponseId,
           color = domain)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  geom_path(alpha = 0.1) +
  geom_smooth(aes(group = capacity), 
              method = "loess", 
              color = "black") +
  scale_x_continuous(breaks = sqrt(as.numeric(levels(factor(d_cap_rating$target_num)))),
                     labels = levels(d_cap_rating$target_ord)) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Perceptions of development, by domain",
       subtitle = "Black line = loess smoothing",
       x = "Target age (numeric, square-root transformed)",
       y = "Response (0 = not at all capable, 100 = completely capable)")

ggplot(d_cap_rating %>%
         arrange(ResponseId, domain, capacity, target_ord) %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_),
           capacity = gsub("_", " ", capacity)),
       aes(x = target_ord, y = response, group = ResponseId,
           color = domain)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  geom_path(alpha = 0.1) +
  geom_smooth(aes(group = capacity), 
              method = "loess", 
              color = "black") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Perceptions of development, by domain",
       subtitle = "Black line = loess smoothing",
       x = "Target age (ordinal)",
       y = "Response (0 = not at all capable, 100 = completely capable)")

Regression models, version 1: Capacity ratings as a function of target age and domain

contrasts_domain_eff_noCOG <- cbind(BOD = c(1, -1, 0, 0),
                                    NEG = c(0, -1, 1, 0),
                                    POS = c(0, -1, 0, 1))
contrasts_domain_eff_noPOS <- cbind(BOD = c(1, 0, 0, -1),
                                    NEG = c(0, 0, 1, -1),
                                    COG = c(0, 1, 0, -1))
contrasts_domain_eff_noNEG <- cbind(BOD = c(1, 0, -1, 0),
                                    POS = c(0, 0, -1, 1),
                                    COG = c(0, 1, -1, 0))
contrasts_domain_eff_noBOD <- cbind(POS = c(-1, 0, 0, 1),
                                    NEG = c(-1, 0, 1, 0),
                                    COG = c(-1, 1, 0, 0))

Target age in months

Effect-coded (comparing all domains to the grand mean, except “bodily sensations”)

Linear effects only

# r1_noBOD <- lmer(response ~ target_num * domain
#                  + (target_num + domain | ResponseId) 
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noBOD))
# saveRDS(r1_noBOD, "./models/r1_noBOD.RDS")
r1_noBOD <- readRDS("./models/r1_noBOD.RDS")
summary(r1_noBOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ target_num * domain + (target_num + domain | ResponseId) +  
    (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 228041.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2123 -0.4827  0.0427  0.5204  4.6348 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept) 195.361  13.977                          
            target_num    9.655   3.107   -0.69                  
            domainPOS    91.876   9.585    0.34 -0.29            
            domainNEG   196.219  14.008    0.57 -0.41 -0.24      
            domainCOG   175.285  13.240   -0.22  0.30 -0.24 -0.58
 capacity   (Intercept) 193.647  13.916                          
            target_num    6.995   2.645   -0.86                  
 Residual               345.942  18.600                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                     Estimate Std. Error t value
(Intercept)           55.2190     5.0015  11.041
target_num             8.1121     0.9584   8.465
domainPOS             -3.8395     8.5475  -0.449
domainNEG              9.0385     8.5719   1.054
domainCOG            -44.3093     8.5670  -5.172
target_num:domainPOS   3.8509     1.6245   2.370
target_num:domainNEG  -2.3475     1.6245  -1.445
target_num:domainCOG   5.4235     1.6245   3.339

Non-linear effects

# r2_noBOD <- lmer(response ~ poly(target_num, 3) * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noBOD))
# saveRDS(r2_noBOD, "./models/r2_noBOD.RDS")
r2_noBOD <- readRDS("./models/r2_noBOD.RDS")
summary(r2_noBOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |  
    ResponseId) + (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 223981

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2929 -0.4858  0.0105  0.5052  5.3522 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept)  89.405   9.455                          
            target_num    9.853   3.139   -0.56                  
            domainCOG   374.649  19.356   -0.47 -0.02            
            domainNEG   435.336  20.865    0.23 -0.48  0.27      
            domainPOS   277.020  16.644    0.04 -0.43  0.50  0.51
 capacity   (Intercept) 193.770  13.920                          
            target_num    7.004   2.646   -0.86                  
 Residual               294.790  17.169                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                               Estimate Std. Error t value
(Intercept)                      66.666      3.914  17.031
poly(target_num, 3)1           2058.398    243.332   8.459
poly(target_num, 3)2           -844.244     17.169 -49.171
poly(target_num, 3)3            407.030     17.169  23.707
domainPOS                         1.595      6.691   0.238
domainNEG                         5.726      6.723   0.852
domainCOG                       -36.656      6.716  -5.458
poly(target_num, 3)1:domainPOS  977.144    412.297   2.370
poly(target_num, 3)2:domainPOS -862.372     29.738 -28.999
poly(target_num, 3)3:domainPOS  507.310     29.738  17.059
poly(target_num, 3)1:domainNEG -595.664    412.297  -1.445
poly(target_num, 3)2:domainNEG  173.963     29.738   5.850
poly(target_num, 3)3:domainNEG  -92.256     29.738  -3.102
poly(target_num, 3)1:domainCOG 1376.179    412.297   3.338
poly(target_num, 3)2:domainCOG   53.207     29.738   1.789
poly(target_num, 3)3:domainCOG -157.616     29.738  -5.300

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Effect-coded (comparing all domains to the grand mean, except “negative emotions”)

Linear effects only

# r1_noNEG <- lmer(response ~ target_num * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noNEG))
# saveRDS(r1_noNEG, "./models/r1_noNEG.RDS")
r1_noNEG <- readRDS("./models/r1_noNEG.RDS")
summary(r1_noNEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ target_num * domain + (target_num + domain | ResponseId) +  
    (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 228041.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2123 -0.4827  0.0427  0.5204  4.6348 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept) 195.361  13.977                          
            target_num    9.655   3.107   -0.69                  
            domainPOS    91.876   9.585    0.34 -0.29            
            domainNEG   196.220  14.008    0.57 -0.41 -0.24      
            domainCOG   175.286  13.240   -0.22  0.30 -0.24 -0.58
 capacity   (Intercept) 193.749  13.919                          
            target_num    6.997   2.645   -0.86                  
 Residual               345.942  18.600                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                     Estimate Std. Error t value
(Intercept)           55.2190     5.0027  11.038
target_num             8.1121     0.9585   8.463
domainBOD             39.1102     8.5568   4.571
domainPOS             -3.8395     8.5497  -0.449
domainCOG            -44.3093     8.5692  -5.171
target_num:domainBOD  -6.9269     1.6248  -4.263
target_num:domainPOS   3.8509     1.6248   2.370
target_num:domainCOG   5.4235     1.6248   3.338

Non-linear effects

# r2_noNEG <- lmer(response ~ poly(target_num, 3) * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noNEG))
# saveRDS(r2_noNEG, "./models/r2_noNEG.RDS")
r2_noNEG <- readRDS("./models/r2_noNEG.RDS")
summary(r2_noNEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |  
    ResponseId) + (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 223981

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2929 -0.4858  0.0105  0.5052  5.3522 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept)  89.406   9.455                          
            target_num    9.853   3.139   -0.56                  
            domainCOG   374.656  19.356   -0.47 -0.02            
            domainNEG   435.338  20.865    0.23 -0.48  0.27      
            domainPOS   277.024  16.644    0.04 -0.43  0.50  0.51
 capacity   (Intercept) 193.741  13.919                          
            target_num    7.003   2.646   -0.86                  
 Residual               294.790  17.169                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       66.666      3.914  17.032
poly(target_num, 3)1            2058.398    243.316   8.460
poly(target_num, 3)2            -844.244     17.169 -49.171
poly(target_num, 3)3             407.030     17.169  23.707
domainBOD                         29.336      6.700   4.378
domainPOS                          1.595      6.691   0.238
domainCOG                        -36.656      6.716  -5.458
poly(target_num, 3)1:domainBOD -1757.660    412.269  -4.263
poly(target_num, 3)2:domainBOD   635.202     29.738  21.360
poly(target_num, 3)3:domainBOD  -257.438     29.738  -8.657
poly(target_num, 3)1:domainPOS   977.144    412.269   2.370
poly(target_num, 3)2:domainPOS  -862.372     29.738 -28.999
poly(target_num, 3)3:domainPOS   507.310     29.738  17.059
poly(target_num, 3)1:domainCOG  1376.179    412.269   3.338
poly(target_num, 3)2:domainCOG    53.207     29.738   1.789
poly(target_num, 3)3:domainCOG  -157.616     29.738  -5.300

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Effect-coded (comparing all domains to the grand mean, except “social abilities & positive emotions”)

Linear effects only

# r1_noPOS <- lmer(response ~ target_num * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noPOS))
# saveRDS(r1_noPOS, "./models/r1_noPOS.RDS")
r1_noPOS <- readRDS("./models/r1_noPOS.RDS")
summary(r1_noPOS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ target_num * domain + (target_num + domain | ResponseId) +  
    (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 228041.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2123 -0.4827  0.0427  0.5204  4.6348 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept) 195.359  13.977                          
            target_num    9.655   3.107   -0.69                  
            domainPOS    91.877   9.585    0.34 -0.29            
            domainNEG   196.219  14.008    0.57 -0.41 -0.24      
            domainCOG   175.285  13.240   -0.22  0.30 -0.24 -0.58
 capacity   (Intercept) 193.641  13.915                          
            target_num    6.994   2.645   -0.86                  
 Residual               345.942  18.600                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                     Estimate Std. Error t value
(Intercept)           55.2190     5.0014  11.041
target_num             8.1121     0.9584   8.465
domainBOD             39.1102     8.5545   4.572
domainNEG              9.0385     8.5718   1.054
domainCOG            -44.3093     8.5669  -5.172
target_num:domainBOD  -6.9269     1.6245  -4.264
target_num:domainNEG  -2.3475     1.6245  -1.445
target_num:domainCOG   5.4235     1.6245   3.339

Non-linear effects

# r2_noPOS <- lmer(response ~ poly(target_num, 3) * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noPOS))
# saveRDS(r2_noPOS, "./models/r2_noPOS.RDS")
r2_noPOS <- readRDS("./models/r2_noPOS.RDS")
summary(r2_noPOS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |  
    ResponseId) + (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 223981

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2929 -0.4858  0.0105  0.5052  5.3522 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept)  89.406   9.455                          
            target_num    9.853   3.139   -0.56                  
            domainCOG   374.652  19.356   -0.47 -0.02            
            domainNEG   435.336  20.865    0.23 -0.48  0.27      
            domainPOS   277.022  16.644    0.04 -0.43  0.50  0.51
 capacity   (Intercept) 193.677  13.917                          
            target_num    7.001   2.646   -0.86                  
 Residual               294.790  17.169                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       66.666      3.914  17.035
poly(target_num, 3)1            2058.398    243.287   8.461
poly(target_num, 3)2            -844.244     17.169 -49.171
poly(target_num, 3)3             407.030     17.169  23.707
domainBOD                         29.336      6.699   4.379
domainNEG                          5.726      6.721   0.852
domainCOG                        -36.656      6.715  -5.459
poly(target_num, 3)1:domainBOD -1757.660    412.216  -4.264
poly(target_num, 3)2:domainBOD   635.202     29.738  21.360
poly(target_num, 3)3:domainBOD  -257.438     29.738  -8.657
poly(target_num, 3)1:domainNEG  -595.664    412.216  -1.445
poly(target_num, 3)2:domainNEG   173.963     29.738   5.850
poly(target_num, 3)3:domainNEG   -92.256     29.738  -3.102
poly(target_num, 3)1:domainCOG  1376.179    412.216   3.338
poly(target_num, 3)2:domainCOG    53.207     29.738   1.789
poly(target_num, 3)3:domainCOG  -157.616     29.738  -5.300

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Effect-coded (comparing all domains to the grand mean, except “cognition & control”)

Linear effects only

# r1_noCOG <- lmer(response ~ target_num * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noCOG))
# saveRDS(r1_noCOG, "./models/r1_noCOG.RDS")
r1_noCOG <- readRDS("./models/r1_noCOG.RDS")
summary(r1_noCOG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ target_num * domain + (target_num + domain | ResponseId) +  
    (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 228041.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2123 -0.4827  0.0427  0.5204  4.6348 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept) 195.360  13.977                          
            target_num    9.655   3.107   -0.69                  
            domainPOS    91.877   9.585    0.34 -0.29            
            domainNEG   196.219  14.008    0.57 -0.41 -0.24      
            domainCOG   175.284  13.239   -0.22  0.30 -0.24 -0.58
 capacity   (Intercept) 193.602  13.914                          
            target_num    6.993   2.644   -0.86                  
 Residual               345.942  18.600                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                     Estimate Std. Error t value
(Intercept)           55.2190     5.0009  11.042
target_num             8.1121     0.9583   8.465
domainBOD             39.1102     8.5536   4.572
domainNEG              9.0385     8.5709   1.055
domainPOS             -3.8395     8.5464  -0.449
target_num:domainBOD  -6.9269     1.6244  -4.264
target_num:domainNEG  -2.3475     1.6244  -1.445
target_num:domainPOS   3.8509     1.6244   2.371

Non-linear effects

# r2_noCOG <- lmer(response ~ poly(target_num, 3) * domain
#                  + (target_num + domain | ResponseId)
#                  + (target_num | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12),
#                  contrasts = list(domain = contrasts_domain_eff_noCOG))
# saveRDS(r2_noCOG, "./models/r2_noCOG.RDS")
r2_noCOG <- readRDS("./models/r2_noCOG.RDS")
summary(r2_noCOG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(target_num, 3) * domain + (target_num + domain |  
    ResponseId) + (target_num | capacity)
   Data: d_cap_rating %>% mutate(target_num = target_num/12)

REML criterion at convergence: 223981

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2929 -0.4858  0.0105  0.5052  5.3522 

Random effects:
 Groups     Name        Variance Std.Dev. Corr                   
 ResponseId (Intercept)  89.405   9.455                          
            target_num    9.853   3.139   -0.56                  
            domainCOG   374.650  19.356   -0.47 -0.02            
            domainNEG   435.333  20.865    0.23 -0.48  0.27      
            domainPOS   277.018  16.644    0.04 -0.43  0.50  0.51
 capacity   (Intercept) 193.787  13.921                          
            target_num    7.004   2.646   -0.86                  
 Residual               294.790  17.169                          
Number of obs: 25896, groups:  ResponseId, 249; capacity, 8

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                       66.666      3.915  17.030
poly(target_num, 3)1            2058.398    243.334   8.459
poly(target_num, 3)2            -844.244     17.169 -49.171
poly(target_num, 3)3             407.030     17.169  23.707
domainBOD                         29.336      6.701   4.378
domainNEG                          5.726      6.723   0.852
domainPOS                          1.595      6.692   0.238
poly(target_num, 3)1:domainBOD -1757.660    412.300  -4.263
poly(target_num, 3)2:domainBOD   635.202     29.738  21.360
poly(target_num, 3)3:domainBOD  -257.438     29.738  -8.657
poly(target_num, 3)1:domainNEG  -595.664    412.300  -1.445
poly(target_num, 3)2:domainNEG   173.963     29.738   5.850
poly(target_num, 3)3:domainNEG   -92.256     29.738  -3.102
poly(target_num, 3)1:domainPOS   977.144    412.300   2.370
poly(target_num, 3)2:domainPOS  -862.372     29.738 -28.999
poly(target_num, 3)3:domainPOS   507.310     29.738  17.059

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Target age in months (square-root transformation)

NOTE: All of these models fail to meet our criterion of random effects being correlation < 0.90. Because of this, I have not included in this print-out and I have not yet modeled non-linear effects after square-root transformation. The next step will be to try simpler random effects models (with fewer random slopes).

Regression models, version 2: Capacity ratings as a function of target age, for each domain separately

NOTE: All of the preregistered models fail to meet our criterion of random effects being correlation < 0.90. Because of this, I instead used simpler random effects models (with fewer random slopes), as indicated in the preregistration.

Target age in months

Bodily sensations

Linear effects only

# r5a_BOD <- lmer(response ~ target_num
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "BOD"))
# saveRDS(r5a_BOD, "./models/r5a_BOD.RDS")
r5a_BOD <- readRDS("./models/r5a_BOD.RDS")
summary(r5a_BOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "BOD")

REML criterion at convergence: 55342

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0547 -0.1326  0.1032  0.2568  3.8043 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 101.283  10.064  
 capacity   (Intercept)   1.223   1.106  
 Residual               106.243  10.307  
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept) 93.31343    1.00055   93.26
target_num   1.32602    0.07675   17.28

Non-linear effects

# r6a_BOD <- lmer(response ~ poly(target_num, 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "BOD"))
# saveRDS(r6a_BOD, "./models/r6a_BOD.RDS")
r6a_BOD <- readRDS("./models/r6a_BOD.RDS")
summary(r6a_BOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "BOD")

REML criterion at convergence: 48865.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1452 -0.2055  0.0208  0.3275  3.8911 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 69.69    8.348   
 capacity   (Intercept)  1.08    1.039   
 Residual               99.44    9.972   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                     Estimate Std. Error t value
(Intercept)            96.002      0.914 105.038
poly(target_num, 3)1  150.369      9.972  15.079
poly(target_num, 3)2 -104.521      9.972 -10.482
poly(target_num, 3)3   74.796      9.972   7.501

Negative emotions

Linear effects only

# r5a_NEG <- lmer(response ~ target_num
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "NEG"))
# saveRDS(r5a_NEG, "./models/r5a_NEG.RDS")
r5a_NEG <- readRDS("./models/r5a_NEG.RDS")
summary(r5a_NEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "NEG")

REML criterion at convergence: 67536.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8357 -0.5649  0.0850  0.7065  3.1869 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 485.5    22.03   
 capacity   (Intercept) 233.9    15.30   
 Residual               573.1    23.94   
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)  63.4263    10.9020   5.818
target_num    5.9054     0.1783  33.128

Non-linear effects

# r6a_NEG <- lmer(response ~ poly(target_num, 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "NEG"))
# saveRDS(r6a_NEG, "./models/r6a_NEG.RDS")
r6a_NEG <- readRDS("./models/r6a_NEG.RDS")
summary(r6a_NEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "NEG")

REML criterion at convergence: 59981.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6632 -0.6339  0.0748  0.6934  3.4866 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 491.7    22.18   
 capacity   (Intercept) 232.9    15.26   
 Residual               549.0    23.43   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                     Estimate Std. Error t value
(Intercept)             72.39      10.89   6.650
poly(target_num, 3)1   731.37      23.43  31.214
poly(target_num, 3)2  -335.14      23.43 -14.304
poly(target_num, 3)3   157.39      23.43   6.717

Social abilities & positive emotions

Linear effects only

# r5a_POS <- lmer(response ~ target_num
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "POS"))
# saveRDS(r5a_POS, "./models/r5a_POS.RDS")
r5a_POS <- readRDS("./models/r5a_POS.RDS")
summary(r5a_POS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "POS")

REML criterion at convergence: 67370.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5267 -0.7058  0.0847  0.7521  2.5604 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 270.5    16.45   
 capacity   (Intercept) 228.1    15.10   
 Residual               572.1    23.92   
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)  50.6737    10.7317   4.722
target_num   11.9782     0.1781  67.256

Non-linear effects

# r6a_POS <- lmer(response ~ poly(target_num, 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "POS"))
# saveRDS(r6a_POS, "./models/r6a_POS.RDS")
r6a_POS <- readRDS("./models/r6a_POS.RDS")
summary(r6a_POS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "POS")

REML criterion at convergence: 58183.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4247 -0.6793 -0.0423  0.6627  3.0011 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 284.4    16.86   
 capacity   (Intercept) 209.4    14.47   
 Residual               420.1    20.50   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                     Estimate Std. Error t value
(Intercept)             68.26      10.29   6.632
poly(target_num, 3)1  1517.77      20.50  74.055
poly(target_num, 3)2  -853.31      20.50 -41.635
poly(target_num, 3)3   457.17      20.50  22.306

Negative emotions

Linear effects only

# r5a_COG <- lmer(response ~ target_num
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "COG"))
# saveRDS(r5a_COG, "./models/r5a_COG.RDS")
r5a_COG <- readRDS("./models/r5a_COG.RDS")
summary(r5a_COG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ target_num + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "COG")

REML criterion at convergence: 62229.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.1523 -0.6120 -0.0335  0.5330  4.1114 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 248.02   15.749  
 capacity   (Intercept)  25.34    5.033  
 Residual               275.14   16.587  
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)  10.7635     3.6912   2.916
target_num   13.4807     0.1235 109.145

Non-linear effects

# r6a_COG <- lmer(response ~ poly(target_num, 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "COG"))
# saveRDS(r6a_COG, "./models/r6a_COG.RDS")
r6a_COG <- readRDS("./models/r6a_COG.RDS")
summary(r6a_COG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ poly(target_num, 3) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "COG")

REML criterion at convergence: 54889.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9686 -0.6187 -0.0333  0.5887  4.1843 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 261.07   16.158  
 capacity   (Intercept)  29.69    5.449  
 Residual               248.48   15.763  
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                     Estimate Std. Error t value
(Intercept)            30.010      3.991   7.519
poly(target_num, 3)1 1717.288     15.763 108.942
poly(target_num, 3)2 -395.519     15.763 -25.091
poly(target_num, 3)3  124.707     15.763   7.911

Target age in months (square-root transformation)

Bodily sensations

Linear effects only

# r7a_BOD <- lmer(response ~ sqrt(target_num)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "BOD"))
# saveRDS(r7a_BOD, "./models/r7a_BOD.RDS")
r7a_BOD <- readRDS("./models/r7a_BOD.RDS")
summary(r7a_BOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "BOD")

REML criterion at convergence: 55196.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.9734 -0.2182  0.0823  0.2948  4.0115 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 101.366  10.068  
 capacity   (Intercept)   1.224   1.106  
 Residual               104.071  10.202  
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       91.6317     1.0088   90.83
sqrt(target_num)   3.6774     0.1733   21.22

Non-linear effects

# r8a_BOD <- lmer(response ~ poly(sqrt(target_num), 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "BOD"))
# saveRDS(r8a_BOD, "./models/r8a_BOD.RDS")
r8a_BOD <- readRDS("./models/r8a_BOD.RDS")
summary(r8a_BOD, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |  
    capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "BOD")

REML criterion at convergence: 48802.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0087 -0.1792 -0.0219  0.3389  4.1311 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 69.73    8.351   
 capacity   (Intercept)  1.08    1.039   
 Residual               98.43    9.921   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                           Estimate Std. Error t value
(Intercept)                 96.0017     0.9141 105.028
poly(sqrt(target_num), 3)1 185.1919     9.9214  18.666
poly(sqrt(target_num), 3)2 -98.7964     9.9214  -9.958
poly(sqrt(target_num), 3)3  36.1618     9.9214   3.645

Negative emotions

Linear effects only

# r7a_NEG <- lmer(response ~ sqrt(target_num)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "NEG"))
# saveRDS(r7a_NEG, "./models/r7a_NEG.RDS")
r7a_NEG <- readRDS("./models/r7a_NEG.RDS")
summary(r7a_NEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "NEG")

REML criterion at convergence: 67281.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8115 -0.5996  0.0806  0.6890  3.5062 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 486.3    22.05   
 capacity   (Intercept) 233.9    15.30   
 Residual               552.6    23.51   
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       57.3409    10.9058   5.258
sqrt(target_num)  14.9239     0.3994  37.370

Non-linear effects

# r8a_NEG <- lmer(response ~ poly(sqrt(target_num), 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "NEG"))
# saveRDS(r8a_NEG, "./models/r8a_NEG.RDS")
r8a_NEG <- readRDS("./models/r8a_NEG.RDS")
summary(r8a_NEG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |  
    capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "NEG")

REML criterion at convergence: 59962.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6297 -0.6396  0.0635  0.6857  3.6271 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 491.8    22.18   
 capacity   (Intercept) 232.9    15.26   
 Residual               547.3    23.39   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                           Estimate Std. Error t value
(Intercept)                   72.39      10.89   6.650
poly(sqrt(target_num), 3)1   809.13      23.39  34.586
poly(sqrt(target_num), 3)2  -164.47      23.39  -7.030
poly(sqrt(target_num), 3)3   -25.42      23.39  -1.087

Social abilities & positive emotions

Linear effects only

# r7a_POS <- lmer(response ~ sqrt(target_num)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "POS"))
# saveRDS(r7a_POS, "./models/r7a_POS.RDS")
r7a_POS <- readRDS("./models/r7a_POS.RDS")
summary(r7a_POS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "POS")

REML criterion at convergence: 65829.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4837 -0.7186  0.0141  0.6928  2.6515 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 274.9    16.58   
 capacity   (Intercept) 228.1    15.10   
 Residual               458.7    21.42   
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       37.4105    10.7351   3.485
sqrt(target_num)  31.2231     0.3638  85.813

Non-linear effects

# r8a_POS <- lmer(response ~ poly(sqrt(target_num), 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "POS"))
# saveRDS(r8a_POS, "./models/r8a_POS.RDS")
r8a_POS <- readRDS("./models/r8a_POS.RDS")
summary(r8a_POS, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |  
    capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "POS")

REML criterion at convergence: 58019

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3200 -0.6911 -0.0292  0.6444  3.3365 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 284.8    16.88   
 capacity   (Intercept) 209.4    14.47   
 Residual               409.1    20.23   
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                           Estimate Std. Error t value
(Intercept)                   68.26      10.29   6.632
poly(sqrt(target_num), 3)1  1736.76      20.23  85.868
poly(sqrt(target_num), 3)2  -540.07      20.23 -26.702
poly(sqrt(target_num), 3)3   -30.24      20.23  -1.495

Cognition & control

Linear effects only

# r7a_COG <- lmer(response ~ sqrt(target_num)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "COG"))
# saveRDS(r7a_COG, "./models/r7a_COG.RDS")
r7a_COG <- readRDS("./models/r7a_COG.RDS")
summary(r7a_COG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ sqrt(target_num) + (1 | ResponseId) + (1 | capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "COG")

REML criterion at convergence: 61697.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.8211 -0.6330 -0.0502  0.5957  3.9643 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 248.79   15.773  
 capacity   (Intercept)  25.34    5.034  
 Residual               254.98   15.968  
Number of obs: 7254, groups:  ResponseId, 279; capacity, 2

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       -0.5581     3.6968  -0.151
sqrt(target_num)  31.4079     0.2713 115.783

Non-linear effects

# r8a_COG <- lmer(response ~ poly(sqrt(target_num), 3)
#                  + (1 | ResponseId)
#                  + (1 | capacity),
#                  d_cap_rating %>%
#                    mutate(target_num = target_num/12) %>%
#                    filter(domain == "COG"))
# saveRDS(r8a_COG, "./models/r8a_COG.RDS")
r8a_COG <- readRDS("./models/r8a_COG.RDS")
summary(r8a_COG, corr = F)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ poly(sqrt(target_num), 3) + (1 | ResponseId) + (1 |  
    capacity)
   Data: 
d_cap_rating %>% mutate(target_num = target_num/12) %>% filter(domain ==  
    "COG")

REML criterion at convergence: 54905.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.8884 -0.6120 -0.0278  0.5740  4.1611 

Random effects:
 Groups     Name        Variance Std.Dev.
 ResponseId (Intercept) 261.04   16.157  
 capacity   (Intercept)  29.69    5.449  
 Residual               249.14   15.784  
Number of obs: 6474, groups:  ResponseId, 249; capacity, 2

Fixed effects:
                           Estimate Std. Error t value
(Intercept)                  30.010      3.991   7.519
poly(sqrt(target_num), 3)1 1756.660     15.784 111.293
poly(sqrt(target_num), 3)2  108.002     15.784   6.842
poly(sqrt(target_num), 3)3 -139.527     15.784  -8.840

Plots: Ratings for developmental factors, by domain and capacity

ggplot(d_dev_factor_rating %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate_at(vars(capacity, dev_factor),
                   funs(gsub("_", " ", .))) %>%
         mutate(dev_factor = factor(
           dev_factor,
           levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
       aes(x = dev_factor, y = response, color = domain)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
  geom_pointrange(data = . %>% group_by(domain, capacity, dev_factor) %>%
                    multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  color = "black", fatten = 4) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Ratings of developmental factors (Version 1)",
       subtitle = "Error bars are bootstrapped 95% CIs",
       x = "Developmental factor", 
       y = "Response (0 = plays no role, 6 = plays a very important role)")

ggplot(d_dev_factor_rating %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate_at(vars(capacity, dev_factor),
                   funs(gsub("_", " ", .))) %>%
         mutate(dev_factor = factor(
           dev_factor,
           levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
       aes(x = reorder(capacity, as.numeric(domain)), 
           y = response, color = domain)) +
  facet_wrap(~ dev_factor, ncol = 5) +
  geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
  geom_pointrange(data = . %>% group_by(domain, capacity, dev_factor) %>%
                    multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  color = "black", fatten = 4) +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1, size = 3))) +
  labs(title = "Ratings of developmental factors (Version 2)",
       subtitle = "Error bars are bootstrapped 95% CIs",
       color = "Domain",
       x = "Capacity (by domain)", 
       y = "Response (0 = plays no role, 6 = plays a very important role)")

ggplot(d_dev_factor_rating %>%
         mutate(dev_factor = factor(
           gsub("_", " ", dev_factor),
           levels = gsub("_", " ", levels(d_dev_factor_rating$dev_factor)))),
       aes(x = dev_factor, y = response)) +
  geom_jitter(alpha = 0.01, height = 0.3, width = 0.3, color = "blue") +
  geom_pointrange(data = . %>% group_by(dev_factor) %>%
                    multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  color = "black", fatten = 2) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Ratings of developmental factors (Version 3)",
       subtitle = "Error bars are bootstrapped 95% CIs",
       x = "Developmental factor", 
       y = "Response (0 = plays no role, 6 = plays a very important role)")

Exploration: EFA/clustering of developmental factors

It could be nice to reduce the number of “devleopmental factors” from 10 down to something more manageable for analysis. Here I explore EFA and hierarchical clustering as dimension reduction methods, and re-plot with these results in mind.

temp <- d_dev_factor_rating %>%
  unite(key, ResponseId, domain, capacity) %>%
  spread(dev_factor, response) %>%
  column_to_rownames("key")

reten_fun(temp, "varimax")

[1] 2

# fa.parallel(temp)
# VSS(temp)

temp_efa <- fa(temp,
               # nfactors = 4,
               nfactors = reten_fun(temp, "varimax"),
               rotate = "varimax")

heatmap_fun(temp_efa) + 
  guides(fill = guide_colorbar(barheight = 8, barwidth = 0.5))

Joining, by = "capacity"
Joining, by = "factor"

library(ggdendro)
library(dendextend)

---------------------
Welcome to dendextend version 1.8.0
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
Or contact: <tal.galili@gmail.com>

    To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------


Attaching package: ‘dendextend’

The following object is masked from ‘package:ggdendro’:

    theme_dendro

The following object is masked from ‘package:stats’:

    cutree

temp_hclust <- temp %>%
  t() %>%
  dist() %>%
  hclust()
temp_hclust_order <- data.frame(order = as.numeric(temp_hclust$order), 
                                dev_factor = as.character(temp_hclust$labels))
temp_hclust %>%
  ggdendrogram(rotate = F) +
  theme_minimal() +
  theme(#axis.title = element_blank(),
        # axis.text.y = element_blank(), 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(title = "Hierarchical agglomerative clustering", 
       subtitle = "Complete linkage (default for stats::hclust() function)",
       x = "Developmental factor", y = "Height")

  # as.dendrogram() %>%
  # set("labels_col", k = 2) %>%
  # plot(horiz = T, xlim = c(-100, 200))

temp2 <- d_dev_factor_rating %>%
  mutate(dev_factor_cluster = case_when(
    dev_factor %in% c("experiments", "people_teach", "brain_changes", 
                      "observes_objects", "observes_people", 
                      "interacts_people") ~ "external",
    dev_factor %in% c("preprogrammed", "womb_experiences", "body_grows", 
                      "senses_improve") ~ "internal",
    TRUE ~ NA_character_)) %>%
  mutate(dev_factor_cluster = factor(dev_factor_cluster),
         response_cent = response - 3) %>%
  left_join(temp_hclust_order)

Joining, by = "dev_factor"
Column `dev_factor` joining factors with different levels, coercing to character vector

contrasts_cluster_eff <- cbind(EXT = c(1, -1))

ggplot(temp2 %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate_at(vars(capacity, dev_factor),
                   funs(gsub("_", " ", .))),
       aes(x = reorder(dev_factor, as.numeric(dev_factor_cluster)), 
           y = response, color = domain, shape = dev_factor_cluster)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  geom_jitter(alpha = 0.025, height = 0.2, width = 0.2, show.legend = F) +
  geom_pointrange(data = . %>% group_by(domain, capacity, 
                                        dev_factor_cluster, dev_factor) %>%
                    multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  color = "black", fatten = 4) +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Ratings of developmental factors (Version 1)",
       subtitle = "Error bars are bootstrapped 95% CIs",
       shape = "Type of developmental factor (per EFA/clustering)",
       x = "Developmental factor", 
       y = "Response (0 = plays no role, 6 = plays a very important role)")

ggplot(temp2 %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate_at(vars(capacity),
                   funs(gsub("_", " ", .))),
       aes(x = reorder(capacity, as.numeric(domain)), 
           y = response, color = domain)) +
  facet_wrap(dev_factor_cluster ~ dev_factor, ncol = 6) +
  geom_jitter(alpha = 0.025, height = 0.2, width = 0.2) +
  geom_pointrange(data = . %>% group_by(domain, capacity, 
                                        dev_factor_cluster, dev_factor) %>%
                    multi_boot_standard(col = "response", na.rm = T),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  color = "black", fatten = 4) +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1, size = 3))) +
  labs(title = "Ratings of developmental factors (Version 3)",
       subtitle = "Error bars are bootstrapped 95% CIs",
       color = "Domain",
       x = "Capacity (by domain)", 
       y = "Response (0 = plays no role, 6 = plays a very important role)")

r_temp <- lmer(response ~ domain * dev_factor_cluster + 
                 # (domain + dev_factor_cluster | ResponseId) +
                 # including both random intercepts yields corr > 0.9
                 (1 + dev_factor_cluster | ResponseId) +
                 (1 | capacity) + (1 | dev_factor),
               temp2,
               contrasts = list(domain = contrasts_domain_eff_noBOD,
                                dev_factor_cluster = contrasts_cluster_eff))

summary(r_temp)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ domain * dev_factor_cluster + (1 + dev_factor_cluster |  
    ResponseId) + (1 | capacity) + (1 | dev_factor)
   Data: temp2

REML criterion at convergence: 87685.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4691 -0.7172  0.0083  0.7176  3.4754 

Random effects:
 Groups     Name                       Variance Std.Dev. Corr 
 ResponseId (Intercept)                0.72087  0.8490        
            dev_factor_clusterinternal 0.31234  0.5589   -0.16
 dev_factor (Intercept)                0.46942  0.6851        
 capacity   (Intercept)                0.08958  0.2993        
 Residual                              2.81132  1.6767        
Number of obs: 22320, groups:  
ResponseId, 279; dev_factor, 10; capacity, 8

Fixed effects:
                                Estimate Std. Error t value
(Intercept)                      3.13801    0.25063  12.520
domainPOS                        0.56263    0.18435   3.052
domainNEG                       -0.22642    0.18435  -1.228
domainCOG                        0.22527    0.18435   1.222
dev_factor_clusterEXT            0.35587    0.22206   1.603
domainPOS:dev_factor_clusterEXT  0.34522    0.01984  17.401
domainNEG:dev_factor_clusterEXT -0.10109    0.01984  -5.095
domainCOG:dev_factor_clusterEXT  0.83940    0.01984  42.309

Correlation of Fixed Effects:
            (Intr) dmnPOS dmnNEG dmnCOG d__EXT dPOS:_ dNEG:_
domainPOS    0.000                                          
domainNEG    0.000 -0.333                                   
domainCOG    0.000 -0.333 -0.333                            
dv_fctr_EXT -0.179  0.000  0.000  0.000                     
dmPOS:__EXT  0.000 -0.022  0.007  0.007  0.000              
dmNEG:__EXT  0.000  0.007 -0.022  0.007  0.000 -0.333       
dmCOG:__EXT  0.000  0.007  0.007 -0.022  0.000 -0.333 -0.333

(Note that domain was effect-coded and leaves out the bodily sensations domain here, but we could re-code it in all the different ways to get the contrast between bodily sensations and the grand mean as desired.)

I see three take-aways from these regression results, all of which I think I can see reflected in the plots above:

1. Participants perceived the domain of social connection and positive emotions to be particularly dramatically shaped by the set of developmental factors we included in this study (relative to the grand mean collapsing across domains).
2. On the whole, participants endorsed “external” and “internal” developmental factors equally strongly. If we were mapping this onto some sort of nature-nuruture distinction (not perfect, but tempting), I think we’d say that there’s no evidence for participants favoring either nature or nurture over the other as an explanation for growth in the capacities we included in this study.
3. However, their endorsement of “external” and “internal” developmental factors varied by domain. They perceived “external” factors to play a larger role in the domains of social connection and positive emotions and especially cognition and control, and a relatively small role in the domain of negative emotions. (Looking at the plots, it seems like external factors were perceived as much more important than internal factors for cognition and control, a bit more important than internal factors for social connection and positive emotions, roughly equally important to internal factors for negative emotions, and less important than internal factors for bodily sensations. These interaction terms just confirm that the difference between “external” and “internal” varied significantly across domains.)

Table: Answers to (optional) free-response question about other developmental factors

d_dev_factor_other %>%
  mutate(domain = recode_factor(
    domain,
    "BOD" = "Bodily sensations",
    "NEG" = "Negative emotions",
    "POS" = "Social abilities & positive emotions",
    "COG" = "Cognition & control",
    .default = NA_character_),
    capacity = gsub("_", " ", capacity)) %>%
  filter(response != "idk") %>%
  distinct(domain, capacity, response) %>%
  arrange(domain, capacity, response) %>%
  kable(caption = "Free responses to optional question about other developmental factors that might play a role (excluding answers like 'N/A,' 'no,' 'not that I can think of', etc.)") %>%
  kable_styling() %>%
  collapse_rows(1:2)

Free responses to optional question about other developmental factors that might play a role (excluding answers like 'N/A,' 'no,' 'not that I can think of', etc.)

domain	capacity	response
Bodily sensations	feeling pain	Again, you're born with this ability
		I think feeling pain is innate
		I think pain is one of the first sensations that humans are able to actually comprehend
		It is just a natural sensation. Our body is programmed to feel bad to avoid further injury or notifiy us that one has occured
		Nervous system
		None come to mind
		Not that I can think of. I think everyone is capable of feeling pain, even at birth, as it is something physical. Whether they understand why they feel it though is a different story
		Pain is likely the very first thing we all experience when we are born. Childbirth is not exactly pleasant for the mother or the baby
		This is another one I would use the word "primal" to describe
		Watching other's reaction to painful situations
		Whether or not they could ignore it or recall feeling pain before and cope
	getting hungry	All children are capable of being hungry
		Biological Need
		child is controlling himself
		everyone gets hungry no matter what the age is . it depends on how we show we are hungry
		evolution, survival mechanisms
		Getting hungry is one of the first things a newborn will experience
		I think we're all just born with the feeling of hunger because he have to have nourishment
		It's just survival at this point
		just something theyre born with
		None come to mind
		the child being active and using enagy
		The child has emotions that cause them to increase their desire to do this
		the child just gets hungry- human nature
		the more experience they have with food, the bigger the capacity for hunger
		We are just preprogrammed to know when we are hungry.
		Yeah, I think this one could be covered with just one answer and needs no further thought or explanation
		You're born with this ability. Technically, you have this ability BEFORE you're born
Negative emotions	feeling distressed	again this is similar to the first one as they age they become more independent and spend less time distressed
		Babies show that they can be distressed from birth, or even during labor.
		Distress is a biological survival mechanism
		Distress is easily felt regardless of age. The factors change. A newborn gets distressed if they are hungry, or tired, or sick, etc. A 5 year old can get distressed if their favorite toy breaks
		Feeling distressed is natural as the child has physical needs and cannot attend to them on their own
		General living
		I believe the feeling of distress is another innate ability that we are able to feel from the moment we're born. I think living, breathing creature that has feelings is the same way from birht
		I do not imagine that any child of any age does not have this ability
		I feel like feeling distressed can be attributed to human's evolution. Distress would be hard coded in the DNA (panic when something bad happens help survival)
		I think babies crying shows distress
		I think it's similar to the pain feeling. They can feel distressed, but they may not understand the feeling
		I think that feeling distressed is somewhat primal. Flight or flight, so to speak. I don't know if I would use the word "biological" to describe it, given my choice, but it's the closest out of the option provided. Then there are layers of socialization, etc. beyond that which amplify it
		i think they always know distress in my opinion
		if the child is left alone or isnt care for
		It's an automatic response
		It's the result of biological function, like pain response
		Maybe hereditary?
		Maybe the amount of nutrition a child gets throughout their life cycle
		None come to mind
		There could be physiological elements that prevent a child from being able to feel distress as acutely as other children
		They are alive
		They can feel stress from the people around them especially parents even inside the womb
		You're born with all emotions, even if you don't understand them yet
	feeling helpless	Again, you're born with emotions. You might not fully comprehend them, but you have the ability to feel them, nontheless
		as they age there is a tendency to become more independent and at that age less aware of danger
		Feeling helpless is something that is learned through experience, not an innate ability
		Feeling helplessness is not really an ability at all, it's an emotion. Most people are capable of feeling helpless
		I don't think you can feel helpless until you become completely self aware, and I'm not sure when that happens
		I feel like a kid would need to have some understanding of their emotions
		I feel that it changes. A younger child may not know or understand the feeling of helpless. An older child my feel they are invinsible or independent.
		I think helplessness stems from a fundamental fear of losing something we've got or not getting something we want (regardless of how well we are able to articulate or understand what that is in a sophisticated way)
		I think it is just mental development. Young children don't really feel self pity
		I think that it all surrounds object perception/realization. Once the child realizes that "food" makes them not hungry or "moving" makes them not hurt
		I think this is something that happens later in life
		My zero is because at the beginning, there's no feeling of self, so you can't feel helpless if you don't feel the self (in my opinion)
		None come to mind
		Not being taken care of or given what it wants would eventually make them real familiar with feeling helpless
		The child has negative experiences that increase this
		The child notices when it is not cared for properly by feelings of hunger,pain etc
		The child's need to be loved and cared for contributes to their feeling helpless, since they are generally moreso helpless in the earlier stages of life than later
		The child's needs change
		This is also hard to rate but I go with my same rationale for happiness. Infants ARE helpless and I think on some leel they know it
Social abilities & positive emotions	feeling happy	A child that feels love feels happiness
		Again, I think happiness (albeit, not amusement) is something relatively innate related to security & safety. I feel like this should be an incredibly easy question to pinpoint to a specific cause, and yet somehow it is not at all
		As the brain grows, we are capable of sensating new emotions.
		I think it's innate in children..Sometimes they interpret it differently as they get older though
		I think its just within them to feel happy
		it is a normal emotion
		None come to mind
		Other external stimuli encourages this
		social environment
		The child becomes more able to notice their own feelings
		the child senses the love shared around him\her
		The child's interests/needs change
		This is a harder one for me to rate. I would say some form of emotional range is present at birth. It would have to be. Beyond that I do not know
		Well fed or not
		yes it depends if the child is in a good home or a bad home
		You're born with this ability
	learning from other people	after the age of 2 or 3 the ability to learn langauges begins to diminish hence my ratings after that age
		I feel like this is just the nature of human beings to learn from the others around (we wouldn't be where we are without it)
		I think there are a ton of factors that lend to learning. I think it was one of the major purposes of the brain and so many, many parts of the brain and body contribute to it
		kids love to copy people they are around and that how they learn
		Learning is an innate ability. Generally, everyone is born with the ability to learn
		None come to mind
		Relationship with parents
		The amount of nutrition the child gets throughout their life cycle
		The more they age, the more connections the brain makes.
Cognition & control	controlling their emotions	A 0 to 9 month old infant is never wrong. You don't tell a 0 to 9 month old infant "No" or "Not right now". After 9 months, then you can start having that conversation
		Children imitate other children
		having very good role models and a good home to learn these in
		I don't think anyone can control thier emotions (make yourself sad, NOW!), only the expression of them
		I know a lot of adults that can't control their emotions. This is something you either have or don't have,
		I would say the type of parenting being done by the parents
		I'm afraid I don't have much to add here. Most of your questions covered it. I would say their environment, but that was, more or less, covered in the questions
		Learning from adults reaction to behavior
		Learning that outcomes can change depending on the the emotion they display
		None come to mind
		overall social environment
		The child has life experiences after birth that increase their ability to do this
		They aren't necessarily learning to control the emotion so much as the expression
		They brains haven't fully developed they cab;'t help themselves and are still exploring their emotions
		They react to positive or negative responses from other humans, if they control, or don't control their emotions
		watching tv or movies
	reasoning about things	Again, up to about 9 months old the wiring is just not there. It beginbs to develop but... to be honest reasoning is something many adults fail at
		Genetic driven intelligence ability
		injury/disability
		It's mostly just natural mental development
		Maybe the amount of food/nutrition they get throughout their life cycle
		None come to mind
		not sure if they can have those skills
		Nuerological development can impact this..for example if they develop brain injuries...autism etc
		The more they age, the more they comprehend and understand.
		watching tv or videos

Plots: Most important developmental factor, by domain and capacity

ggplot(d_dev_factor_most_important_choice %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate(dev_factor = case_when(
           grepl("Other", response) ~ "other",
           grepl("teach", response) ~ "people teach",
           grepl("experiments", response) ~ "experiments",
           grepl("womb", response) ~ "womb experiences",
           grepl("interacts", response) ~ "interacts people",
           grepl("preprogrammed", response) ~ "preprogrammed",
           grepl("objects", response) ~ "observes objects",
           grepl("observes the people", response) ~ "observes people",
           grepl("body grows", response) ~ "body grows",
           grepl("brain changes", response) ~ "brain changes",
           grepl("senses improve", response) ~ "senses improve"),
           dev_factor = factor(
             dev_factor,
             levels = c(gsub("_", " ", 
                             levels(d_dev_factor_rating$dev_factor)),
                        "other"))) %>%
         mutate_at(vars(capacity, response),
                   funs(gsub("_", " ", .))),
       aes(x = dev_factor, fill = domain)) +
  facet_wrap(~ domain ~ capacity, ncol = 4) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  geom_bar() +
  labs(title = "Choice of most important developmental factor (Version 1)",
       x = "Developmental factor", 
       y = "Number of participants")

temp_other <- d_dev_factor_most_important_choice %>% 
  filter(grepl("Other", response)) %>% 
  count(domain, capacity)
ggplot(d_dev_factor_most_important_choice %>%
         mutate(domain = recode_factor(
           domain,
           "BOD" = "Bodily sensations",
           "NEG" = "Negative emotions",
           "POS" = "Social abilities & positive emotions",
           "COG" = "Cognition & control",
           .default = NA_character_)) %>%
         mutate(dev_factor = case_when(
           grepl("Other", response) ~ "other",
           grepl("teach", response) ~ "people teach",
           grepl("experiments", response) ~ "experiments",
           grepl("womb", response) ~ "womb experiences",
           grepl("interacts", response) ~ "interacts people",
           grepl("preprogrammed", response) ~ "preprogrammed",
           grepl("objects", response) ~ "observes objects",
           grepl("observes the people", response) ~ "observes people",
           grepl("body grows", response) ~ "body grows",
           grepl("brain changes", response) ~ "brain changes",
           grepl("senses improve", response) ~ "senses improve"),
           dev_factor = factor(
             dev_factor,
             levels = c(gsub("_", " ", 
                             levels(d_dev_factor_rating$dev_factor)),
                        "other"))) %>%
         mutate_at(vars(capacity, response),
                   funs(gsub("_", " ", .))) %>%
         filter(dev_factor != "other"),
       aes(x = reorder(capacity, as.numeric(domain)), fill = domain)) +
  facet_wrap(~ dev_factor, ncol = 5) +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  geom_bar() +
  labs(title = "Choice of most important developmental factor (Version 2)",
       subtitle = paste0("Excluding answers of 'other' (n < ", 
                         max(temp_other$n) + 1, ") for any capacity"),
       color = "Domain",
       x = "Capacity (by domain)", 
       y = "Number of participants")

d_dev_factor_most_important_free %>%
  mutate(domain = recode_factor(
    domain,
    "BOD" = "Bodily sensations",
    "NEG" = "Negative emotions",
    "POS" = "Social abilities & positive emotions",
    "COG" = "Cognition & control",
    .default = NA_character_),
    capacity = gsub("_", " ", capacity)) %>%
  distinct(domain, capacity, response) %>%
  arrange(domain, capacity, response) %>%
  kable(caption = "Free responses when 'other' was selected as most important developmental factor") %>%
  kable_styling() %>%
  collapse_rows(1:2)

Free responses when 'other' was selected as most important developmental factor

domain	capacity	response
Bodily sensations	feeling pain	everyone is capable of feeling pain.
	getting hungry	All children are capable of being hungry.
Negative emotions	feeling distressed	As soon as his consious develops they are considered alive and capable of feeling.
		emotional development
		most id not all children are capable of experinceing distress not matter what.
	feeling helpless	Experiences in life
		I really don't know.
		Not sure helplessness is felt so young
Social abilities & positive emotions	feeling happy	Every child is capable of feeling happiness no matter what.
		it is just a normal response too something enjoyed
	learning from other people	not sure
Cognition & control	controlling their emotions	impossibe

Demographics

ggplot(d_demo, aes(x = Duration/60)) + 
  geom_histogram(binwidth = 2) +
  geom_vline(xintercept = median(d_demo$Duration/60), color = "blue", lty = 2) +
  scale_x_continuous(breaks = seq(0, 10000, 4)) +
  labs(title = "Duration of study (according to Qualtrics)",
       subtitle = "Blue dotted line marks median",
       x = "Duration (in minutes)",
       y = "Number of participants")

ggplot(d_demo, aes(x = Age)) + 
  geom_histogram(binwidth = 2) +
  geom_vline(xintercept = median(d_demo$Age), color = "blue", lty = 2) +
  scale_x_continuous(breaks = seq(0, 10000, 4)) +
  labs(title = "Particpiant age (self-reported)",
       subtitle = "Blue dotted line marks median",
       x = "Age (in years)",
       y = "Number of participants")

ggplot(d_demo, aes(x = GenderSex)) + 
  geom_bar() +
  labs(title = "Particpiant gender/sex (self-reported)",
       x = "Gender/sex",
       y = "Number of participants")

ggplot(d_demo, aes(x = gsub('(.{1,30})(\\s|$)', '\\1\n', RaceEthnicity_collapse))) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Particpiant race/ethnicity (self-reported)",
       x = "Race/ethnicity",
       y = "Number of participants")

ggplot(d_demo, aes(x = FirstLang)) + 
  geom_bar() +
  labs(title = "Particpiant first language (self-reported)",
       x = "Language",
       y = "Number of participants")

ggplot(d_demo, aes(x = factor(Education,
                              levels = levels(d$Education),
                              labels = gsub('(.{1,30})(\\s|$)', '\\1\n', 
                                            levels(d$Education))))) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Particpiant educational attainment (self-reported)",
       x = "Highest level of education completed",
       y = "Number of participants")

ggplot(d_demo, aes(x = Income)) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Particpiant household income (self-reported)",
       x = "Annual household income",
       y = "Number of participants")

ggplot(d_demo, aes(x = HouseholdSize)) + 
  geom_histogram(binwidth = 1) +
  geom_vline(xintercept = median(d_demo$HouseholdSize), color = "blue", lty = 2) +
  scale_x_continuous(breaks = seq(0, 10000, 1)) +
  labs(title = "Particpiant household size (self-reported)",
       subtitle = "Blue dotted line marks median",
       x = "Number of people in household (adults and children)",
       y = "Number of participants")

ggplot(d_demo, aes(x = MaritalStatus)) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Particpiant marital status (self-reported)",
       x = "Marital status",
       y = "Number of participants")

ggplot(d_demo, aes(x = Parent)) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Particpiant parent status (self-reported)",
       subtitle = "'NA' indicates response of 'Prefer not to say'",
       x = "Parent status",
       y = "Number of participants")

ggplot(d_demo %>% filter(Parent == "Yes"), aes(x = ChildrenNumber)) + 
  geom_histogram(binwidth = 1) +
  geom_vline(xintercept = median(d_demo[d_demo$Parent == "Yes",]$ChildrenNumber, na.rm = T), 
             color = "blue", lty = 2) +
  scale_x_continuous(breaks = seq(0, 10000, 1)) +
  labs(title = "Number of children among parents (self-reported)",
       subtitle = "Blue dotted line marks median",
       x = "Number of children (among parents)",
       y = "Number of participants")

ggplot(d_demo %>% filter(Parent == "Yes"), 
       aes(x = factor(ChildrenYoungestAge_collapse,
                      levels = levels(d_demo$ChildrenYoungestAge_collapse),
                      labels = gsub('(.{1,30})(\\s|$)', '\\1\n', 
                                    levels(d_demo$ChildrenYoungestAge_collapse))))) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Age of youngest child among parents (self-reported)",
       x = "Age of child",
       y = "Number of participants")

ggplot(d_demo %>% filter(Parent == "Yes"), 
       aes(x = factor(ChildrenOldestAge_collapse,
                      levels = levels(d_demo$ChildrenOldestAge_collapse),
                      labels = gsub('(.{1,30})(\\s|$)', '\\1\n', 
                                    levels(d_demo$ChildrenOldestAge_collapse))))) + 
  geom_bar() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Age of oldest child among parents (self-reported)",
       x = "Age of child",
       y = "Number of participants")